Demonstration of the new CLiMAF operator 'plotmap',

See doc at https://climaf.readthedocs.io/en/master/scripts/plotmap.html

plotmap is dedicated to replace operator 'plot'; it is backed by script plotmap.py which uses matplotlib, cartopy and geocat

from climaf.api import *

python => 3.10.13 | packaged by conda-forge | (main, Oct 26 2023, 18:07:37) [GCC 12.3.0]
---
Required softwares to run CliMAF => you are using the following versions/installations:

CliMAF climaf_version = 3.0
CliMAF install => /home/ssenesi/climaf_installs/climaf_running
Cache directory set to : /data/ssenesi/climaf_cache (use $CLIMAF_CACHE if set) 
Cache directory for remote data set to : /data/ssenesi/climaf_cache/remote_data (use $CLIMAF_REMOTE_CACHE if set) 
warning  : defining curl_tau_atm : command ferret is not executable
warning  : Binary cdftools not found. Some operators won't work
Available macros read from ~/.climaf.macros are : []
error    : second argument ('curl_tau_atm') must be a script or operator, already declared

# Some CMIP data
data = fds(
    cpath + "/../examples/data/tas_Amon_CNRM-CM5_historical_r1i1p1_185001-185212.nc")

# More example data
dg = ds(project="example", simulation="AMIPV6ALB2G",
        variable="tas", period="1981", frequency="monthly")

# When developing the script of a CliMAF operator, we may wish to systematically
# recompute the object we look at
# Redefining function cshow allows for that
def cshow(obj):
    climaf.api.cshow(obj, drop=True)

Basics

# For reference, a figure using the old script (gplot.ncl)
cshow(plot(data, title="Reference, using gplot.ncl"))

# Same with new script. Mapping values to colors is slightly different 
cshow(plotmap(data, title="Using plotmap"))

# One can easily tune the number of levels used
cshow(plotmap(data, title="hand-tuned levels #", levels=10))

# Changing the Reference longitude
cshow(plotmap(data, title="Changing the Reference longitude", 
              proj="PlateCarree", proj_options={"central_longitude":90}))

# The part of the coordinates range shown is called the 'extent'
# You can change it, but you must specify longitudes that fits in central_longitude +/- 180°
# So, by default (i.e. central_longitude=180), that fits in  [ 0, 360 [
cshow(plotmap(data,axis_methods={'set_extent': {'extents': (0 , 90, -10, 30) }}))

Title and other strings are managed by function geocat.viz.util.set_titles_and_labels wich uses arguments in 'title_options'

# Playing with title size
cshow(plotmap(data, title='Main title',
              title_options = {'maintitlefontsize':50}))

# If you provide both 'main_title' in 'title_options' and 'title', the former prevails
cshow(plotmap(data, title='Main title',
              title_options = {'maintitle':'title_options prevails'}))

# Controling upper left and right strings
cshow(plotmap(data, title='Main title',
              title_options = {
                  'lefttitle':'Left String',
                  'righttitle':'Right String',
                  'righttitlefontsize': 6,
              }))

# Alas, using a subtitle ruins the vertical spacing (this is a bug in a geocat-viz function)
cshow(plotmap(data, title='Main title',
              title_options = {
                  'lefttitle':'Left String',
                  'righttitle':'Right String',
                  'righttitlefontsize': 6,
                  'subtitle':'Subtitle',}))

Playing with colors

# Play with number of colors
cshow(plotmap(data, levels=15))

# Choosing each color
mycolors = 'white,black,RoyalBlue,LightSkyBlue,PowderBlue,lightseagreen,PaleGreen,Wheat,Brown,Pink'
cshow(plotmap(data, color= mycolors))

# Choosing the mapping of values to colors
mycolors = 'white,black,RoyalBlue,LightSkyBlue,PowderBlue'
mylevels = '210,220,240,300,310'
cshow(plotmap(data, color= mycolors, levels=mylevels))

# 'colors' and levels" can also be python lists
mycolors = [ 'white' , 'black' , 'RoyalBlue' , 'LightSkyBlue' , 'PowderBlue' ]
mylevels = [210,220,240,300,310]
cshow(plotmap(data, color= mycolors, levels=mylevels))

# Adding contours between colored regions
mylevels = [210,220,240,290,300,310]
cshow(plotmap(data, levels=mylevels, contours=1))

# NCL colormaps are available
cshow(plotmap(data, color = 'helix'))

# How to change gridlines look
cshow(plotmap(data, axis_methods={
    'gridlines': { 'linewidths': 5, 'linestyles': 'dashed', 'draw_labels' : {"bottom": "x", "left": "y"}}}))

# To avoid smoothing, use colormap engine 'pcolormesh' 
cshow(plotmap(data, clre="pcolormesh"))

Data selection

# Select on time index
cshow(plotmap(data, time=0, print_time=True))

# Again selection on time index, with different index value
cshow(plotmap(data, time=1, print_time=True))

# Selection on date
cshow(plotmap(data, date='185003', print_time=True))

# Selection by providing explicit 'selection options' arguments
cshow(plotmap(data, clrso={"isel" : {"time": 0}}, print_time=True))

# Again: selection by providing explicit 'selection options' argument here with a date
cshow(plotmap(data, clrso={'sel' : {'time': '1850-02'}}, print_time=True))

Mimicking various gplot options/args

# Changing image size
cshow(plotmap(data, title="Changing image size", resolution="600x400"))

# Showing only a data range. Here is gplot.ncl reference
cshow(plot(data, min=210, max=250, delta=5, title="A range, with 'plot'"))

# Showing only a data range. Here is potmap version
cshow(plotmap(data, min = 210, max = 250, delta=5, title="A range, with 'plotmap'"))

# What if surrounding white space is left
cshow(plotmap(data,trim=False))

# Request an horizontal colorbar
cshow(plotmap(data,vcb=False))

# Focus plot on ocean
cshow(plotmap(data,focus='ocean'))

# Focus plot on land
cshow(plotmap(data,focus='land'))

# Apply offset and change displayed units accordingly
cshow(plotmap(data, offset=-273.15, units="C"))

# Polar stereo projection 
cshow(plotmap(data,proj='NH70'))

# Polylines
cshow(plotmap(data,xpolyline="0 90",ypolyline="0 45", polyline_options={'color': 'blue'}))

# Nemo data :
nd = fds(cpath + "/../examples/data/tos_Omon_CNRM_gn_185001-185003.nc")
cshow(plotmap(nd))

# Filling the gaps near ORCA grid poles needs another plot engine  :
cshow(plotmap(nd, clre="pcolormesh"))

Pattern shading (hatching)

See also hatch style reference

# Pattern-shading up to two fields, provided as 5th and/or 6th plotmap arg
# shdl (resp. shdl2) provides levels between which shading occurs
# shdh (resp shdh2) provides patterns list, which is reused circularly if needed
#    None means no pattern, repeated pattern means increased density
cshow(plotmap(dg, "", "", "", dg, shdl=[210, 240, 270, 285, 300],shdh = ['+', None, '///']))

Vectors

# Vector plots. Get data
ua = fds(cpath + "/../examples/data/uas_CNRM-CM6_sample.nc")
va = fds(cpath + "/../examples/data/vas_CNRM-CM6_sample.nc")

# Default vector representation type is 'quiver' (arrows)
# Vector components are args 3 and 4 of plotmap
cshow(plotmap("","",ua,va))

Tuning vector grid size and arrows attributes. See quiver doc for more options (but don't play with args X,Y, u, v)

cshow(plotmap("","",ua,va,
              vecg=50, 
              veco={'color':'blue', 'headwidth':2.5, 'headlength':4}))

# Controling the map projection
cshow(plotmap("","",ua,va, proj="Mollweide"))

Another vector representation type is barbs. See its doc for more options

cshow(plotmap("","",ua,va,vecty="barbs", vecg=40, veco={'length':3.5, 'barbcolor' : 'blue'},))

The third vector representation type is streamlines. Here is the doc for more options

cshow(plotmap("","",ua,va,vecty="streamplot", veco={'density' : 3 , 'linewidth': 0.7}))

# Superimposing colored map and streamplot
cshow(plotmap(dg,"",ua,va,
              vecty="streamplot",  
              veco={'density' : 3 , 'linewidth': 0.7},
              axis_methods={'set_extent': {'extents': (0. , 60, 0., 60.) }}))

Output formats

# Output in pdf format
f=cfile(plotmap(dg, format='pdf'))
print(f)
! display $f

/data/ssenesi/climaf_cache/d0/1e770f3522e3af4064b534f2bbc505fad6cbf729f5d78a22e7c885.pdf
/bin/bash: /net/nfs/tools/Users/SU/jservon/spirit-2021.11_envs/20240209/lib/libtinfo.so.6: no version information available (required by /bin/bash)

pdftk Warning: unexpected case 1 in LoadDataFile(); continuing

# Output in eps format
f=cfile(plotmap(dg, format='eps'))
print(f)
! display $f

/data/ssenesi/climaf_cache/62/968cb8e251f72c53fa54c57622de27180197b68f2fc06a483442e2.eps
/bin/bash: /net/nfs/tools/Users/SU/jservon/spirit-2021.11_envs/20240209/lib/libtinfo.so.6: no version information available (required by /bin/bash)

# Launching plotmap with format='show' uses 
# matplotlib.show(block=True), which pops up a window. 
# This works only outisde a Notebook 

#cfile(plotmap(dg, format="show"))  

Playing with projections, and using data on a rectangular projected grid

# Next data was generated on a grid that uses Lambert2 projection
l2 = fds(cpath + "/../examples/data/sfcWind_aladin_ext.nc")

# As a default, target projection is PlateCarree
# Plotmap scrutinizes file metadata for identifying data projection parameters
cshow(plotmap(l2))

# We can also explicitly specify data grid's projection parameters
# This is useful if data file metadata is missing or not 'understandable' by plotmap
transform = { 
    'clrt' : 'LambertConformal', 
    'clrto' : {
        'central_longitude' : 10,
        'central_latitude'  : 37,
        'standard_parallels' : (37,37),
    }
    }
cshow(plotmap(l2, **transform))

# Another way to explicitly specify data grid's projection parameters is 
# to provide a file which includes correct metadata
cshow(plotmap(l2, clrt=cfile(l2)))

# This plot engine option speeds up computaton, but may damage the plot
cshow(plotmap(l2, clreo={'transform_first':True}))

# We can request another target projection
cshow(plotmap(l2, proj="Stereographic"))

# If we want to plot on the data grid, we may describe it explicitly, 
# Here we use args 'proj' and 'proj_options', and pack them in a small dict
projection = { 
    'proj' : 'LambertConformal', 
    'proj_options' : {
        'central_longitude' : 10,
        'central_latitude'  : 37,
        'standard_parallels' : (37,37),
    }
    }
# Note : data will actually be remapped on its own grid (To be checked)
cshow(plotmap(l2, **projection))

# But we can also just refer to the projection of the colored field 
cshow(plotmap(l2, proj="clr"))

# We can also describe target projection using a file with relevant metadata
cshow(plotmap(l2, proj=cfile(l2)))

# We can explictly request that the data is not remapped. But there is no watch dog here !
cshow(plotmap(l2, proj=cfile(l2), clrt="no_remap"))

# Still without remapping, and using another plot engine to 'see' data grid cells
cshow(plotmap(l2, proj=cfile(l2), clrt="no_remap", clre="pcolormesh"))

# Selecting the plot domain, using 'extents' (here projection PlateCarree is implicit)
cshow(plotmap(l2, axis_methods={'set_extent': {'extents': (0 , 70, 20, 60) }}))

Advanced features

The methods of cartopy.mpl.geoaxes.GeoAxes can be called using arg axis_methods

cshow(plotmap(dg, axis_methods={'coastlines': {'color': 'white'}}))

cshow(plotmap(dg, axis_methods={'gridlines': {}}))

# Add feature (implictly from cartopy.feature) 
cshow(plotmap(dg, axis_methods={
    'add_feature': {'feature': 'LAND', 'facecolor': 'black', 'zorder': 1}}))

# Add feature on contoured field 
cshow(plotmap("",dg, axis_methods={
    'add_feature': {'feature': 'LAND', 'facecolor': 'black', 'zorder': 1}}))

Also, Geocat-viz methods can be invoked using arg 'gv_method'

cshow(plotmap(dg, 
              gv_methods={'add_major_minor_ticks': { 
                              'labelsize':'small',
                              'x_minor_per_major':2
                              }}))

Some methods of matplotlib.pyplot can be called using arg plt_methods

cshow(plotmap(dg, plt_methods={
    'text': {'x':-120, 'y': 45, 's': 'mytext', 'horizontalalignment': 'left'}}))